echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 9.S (Runge-Kutta Method for a System).
% Section	9.7,	Runge-Kutta Methods, Page 475
echo on; clc; format long; hold off; clear

% This program implements the Runge-Kutta method
% for solving the initial value problem

%     Z' = F(t,Z)     with    Z(a) = Za

% where the vector notation is equivalent to:

%     x' = f(t,x,y)   with    x(a) = xa
%     y' = g(t,x,y)           y(a) = ya

% Formulas for f(t,x,y) and g(t,x,y) are used to form
%	the vector function F(t,Y)	which is stored in  Fn.m 

% function W = fn(t,Z)
% x = Z(1);  y = Z(2);
% W = [(x - x*y - x^2/10), (x*y - y - y^2/20)];

delete fn.m
diary fn.m; disp('function W = fn(t,Z)');...
            disp('x = Z(1);  y = Z(2);');...
            disp('W = [(x - x*y - x^2/10), (x*y - y - y^2/20)];');...
diary off;

fn(0,[0 0]); % Remark. fn.m and rks4 are used in Algorithm 9.S
pause	       % Press any key to continue.

clc;
%	Place the endpoints of [a,b] in  a  and  b.

%	Place the initial value Za = Z(a) in Za.

%	Place the number of subintervals in  m.

a  = 0;
b  = 15;
Za = [2  1];
m  = 150;
[T,Z] = rks4('fn',a,b,Za,m);
P = [T;Z']';
points = P(1:10:length(P),:);

pause	% Press any key to see the list of solution points.

clc;, clg;
X = Z(:,1);
Y = Z(:,2);
a = min(X)-0.3;
b = max(X)+0.3;
c = min(Y)-0.3;
d = max(Y)+0.3;
axis([0 2.1 0 1.8]);...
plot(X,Y,'g');...
hold on;...
plot([0 2.1],[0 0],'b',[0 0],[0 1.8],'b');...
xlabel('X');...
ylabel('Y');...
title('Runge-Kutta solution to Z` = F(t,Z)');...
grid;...
axis;...
hold off;...
shg; pause	% Press any key to continue.

Mx1 = 'Runge-Kutta solution to Z` = F(t,Z).';
Mx2 = '     t(k)               x(k)               y(k)';
clc,echo off,diary output,...
disp(''),disp(Mx1),...
disp(''),disp(Mx2),disp(points),diary off,echo on
